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Abstract 

Recent advance in nanotechnology has led to rapid advances in nanofluidics, which has been established as a 
reliable means for a wide variety of applications, including molecular separation, detection, crystallization and 
biosynthesis. Although atomic and molecular level consideration is a key ingredient in experimental design and 
fabrication of nanfluidic systems, atomic and molecular modeling of nanofluidics is rare and most simulations 
at nanoscale are restricted to one- or two-dimensions in the literature, to our best knowledge. The present work 
introduces atomic scale design and three-dimensional (3D) simulation of ionic diffusive nanofluidic systems. We 
propose a variational multiscale framework to represent the nanochannel in discrete atomic and/or molecular 
detail while describe the ionic solution by continuum. Apart from the major electrostatic and entropic effects, 
the non-electrostatic interactions between the channel and solution, and among solvent molecules are accounted 
in our modeling. We derive generalized Poisson-Nernst-Planck (PNP) equations for nanofluidic systems. 
Mathematical algorithms, such as Dirichlet to Neumann mapping and the matched interface and boundary 
(MIB) methods are developed to rigorously solve the aforementioned equations to the second-order accuracy 
in 3D realistic settings. Three ionic diffusive nanofluidic systems, including a negatively charged nanochannel, 
a bipolar nanochannel and a double-well nanochannel are designed to investigate the impact of atomic charges 
to channel current, density distribution and electrostatic potential. Numerical findings, such as gating, ion 
depletion and inversion, are in good agreements with those from experimental measurements and numerical 
simulations in the literature. 
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channel. 
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1. Introduction 


Nanofluidics refers to the study of the transport of ions and/or molecules in confined solutions as well 
as fluid flow through or past structures with one or more characteristic nanometer dimensions [26] [75]. The 
dramatic advances in microfluidics in the 1990s and the introduction of nanoscience, nanotechnology and 
atomic fabrication in recent years have given its own name to nanofluidics [32] . Nanofluidic systems have been 
extensively exploited for molecule separation and detection, nanosensing, elucidation of complex fluid behavior 
and for the discovery of new physical phenomena that are not observed or less influential in macrofluidic or 
microfluidic systems [82 [. Some of such phenomena include double-layer overlap, ion permitivity, diffusion, 
ion-current rectification, surface charge effect and entropic forces [75l [107], 

One major feature of a nanofluidic system is its structural characteristic. Nanofluidic structures can be 
classified into nanopores and nanochannels and, in fact, these two terms are exchangeable in many cases [107] , 
A nanopore has comparatively short length formed perpendicularly through various materials, such as a bipore 
consisting of proteins, i.e., a -hemolysin and a solid-state pore [SUEZ]. An example of solid-state pore is a 
set of nanopores in a silicon nitride membrane which enables the detection of folding behaviors of a single 
double-stranded DNA [62]. On the other hand, a nanochannel has relatively larger dimensions of depth and 
width, usually fabricated in a planar format, and is often equipped with other sophisticated devices to control 
or influence the transport inside the channel [ 107 . For instance, Perry et al. demonstrates the rectifying 
effect of a funnel-shape nanochannel based on different movements of counterions at its tip and base m a 
nano-scaled channel usually has either a cylindrical or a conical geometry [Ml. In a cylindrical channel, the 
flow direction does not influence on current, but surface charges and applied external voltage alter the flux 
of ions with opposite sign charges. However, the difference in the size of pores in a conical channel brings 
different ionic conductance patterns depending on the flow direction. 

The other major feature of a nanofluidic system is its interactions. It is the interaction at nanoscale that 
distinguishes a nanofluidic system from an ordinary fluid system. Certainly, most interactions are directly 
inherited from the chemical and physical properties of the nanostructure, such as the geometric confinement, 
steric effect, polarization and charge. Some other interactions are controlled by flow conditions, i.e., ion 
composition and concentration, and applied external fields. Therefore, the interactions of a nanofluidic system 
is determined by its structure and flow conditions. The function of a nanofluidic system is in turn determined 
by all the interactions. Usually, most nanofluidic systems do not involve any chemical reactions. In this 
case, steric effects, van der Waals interactions and electrostatic interactions are pivoting factors. Therefore, 
in nanofluidic systems, microscopic interactions dominate the flow behavior, while in macroscopic flows and 
some microfluidics, continuum fluid mechanics governs and microscopic effects are often negligible. Typically, 
microscopic and macroscopic behaviors co-exist in a microfluidic system. 

Characteristic length scales, such as Reynolds number, Biot number and Nusselt number, are important to 
the macroscopic fluid flows. For most nanofluidic systems, one of most important characteristic length scales 
is the Debye length X D = x ^ £ °^ bT 2 , where 6 is the dielectric constant of the solvent, £q is the permittivity 

V ca Q a 

of vacuum, ks is the Boltzmann constant, T is the absolute temperature, and C a o and q a are, respectively, 
the bulk ion concentration and the charge of ion species a [26 . The Debye length describes the thickness (or, 
precisely, ^th reduction) of electrical double layer (EDL). Essentially, ionic fluid behaves like a microscopic 
flow within the EDL region, while acts as a macroscopic flow far beyond the Debye length. By the Gouy- 
Chapman-Stern model, the EDL is divided into three parts: the inner Helmholtz plane, outer Helmholtz 
plane and diffuse layer m • While the inner Helmholtz plane consists of non-hydrated coions and counterions 
that are attached to the channel surface, the outer Helmholtz plane contains hydrated or partially hydrated 
counterions. Moreover, the part between the inner and outer Helmholtz planes is called the Stern layer. 

Note that the EDL applies not only to the layer near the nanochannel, but also the layer around a charged 
biomolecule in the flow. Consequently, many microfluidic devices with quite large channel dimensions exhibit 
microscopic flow characteristic when the fluid consists of large macromolecules and solvent. The possible 
deformation, aggregation, folding and unfolding of the macromolecules in the fluidic system make the fluid flow 
behavior complex and intriguing [40: . Nevertheless, for rigid macromolecules, the effective channel dimensions 
can be estimated by subtracting the macromolecular dimension from the physical dimension of the channel. 
The resulting system may be approximated by simple ions for most analysis. 

Specifically, the charge on the wall surface derives electrostatic interactions and electrokinetic effects when 
ions in a solution are sufficiently close to channel wall |83, [103]. Since the surface-to-volume ratio is excep- 
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tionally high in a nanoscale channel, surface charges induce a unique electrostatic screening region, i.e., EDL 
m- In fact, it attracts ions charged oppositely (counterions) and repels ions having the same charge (coions) 
to sustain the electroneutrality of an aqueous solution confined in a channel. Physically, the EDL region only 
contains bound or mobile counterions and typically covers the nano-sized pore of a channel. Therefore, the 
oppositely charged ions mainly constitute the electrical current through a micro- or nano-channel [75 L f88]. 

The rectification of ionic current, which is one of the distinct transport properties of nanofluidic channels, 
can further elucidate the flow pattern and formation of the fluid through a nanochannel [56]. This phenomenon 
usually occurs when surface charge distribution, applied electric field, bulk concentration and/or channel 
geometry are properly manipulated along the channel axis [21]. Pu et al. conducted experiments to present 
ion-enrichment and ion-depletion effects on nanochannels to show that the rectification begins with these two 
effects m- In their design, an applied field gave rise to accumulation of all ions at the cathode and absence 
of all ions at the anode of the channels. Ion selectivity is another important feature which enables nano-sized 
channels to work as an ionic filter [88]. It is defined as the ratio of the difference between currents of cations and 
anions to the total current delivered by both ions. Vlassiouk and his colleagues examined the ion selectivity 
of single nanometer channels under various conditions including channel dimension, buffer concentration and 
applied voltage [88] . 

Nanofluidics has been extensively studied in chemistry, physics, biology, material science, and many areas 
of engineering m- The primary purpose of most studies is to separate and/or detect biological substances 
in a complex solution [69]. A variety of nanofluidic devices have been produced using extraordinary trans¬ 
port behaviors caused by steric restriction, polarization and electrokinetic principles mm- For instance, a 
nanofluidic diode is an outstanding tool to take the advantage of the rectifying effect of ionic current through 
a nanochannel [56]. The nanofluidic diodes have been developed to govern the flow inside the channel by 
breaking the symmetry in channel geometry, surface charge arrangement and bulk concentration under the 
influence of applied voltage nm]. Additionally, the design and fabrication of nanofluidics for molecular bi¬ 
ology applications is a new interdisciplinary field that makes use of precise control and manipulation of fluids 
at submicrometer and nanometer scales to study the behavior of molecular and biological systems. Because of 
the microscopic interactions, fluids confined at the nanometer scale can exhibit physical behaviors which are 
not observed or insignificant in larger scales. When the characteristic length scale of the fluid coincides with 
the length scale of the biomolecule and the scale of the Debye length, nanofluidic devices can be employed for 
a variety of interesting basic measurements such as molecular diffusion coefficients !H!, enzyme reaction rates 
nmm], pH values [65] [97], and chemical binding affinities [54]. Micro- and nanofluidic techniques have been 
instrumented for polymerase chain reaction (PCR) amplifications [8], macromolecule accumulator [22] 198] . 
electrokinetics gnnn], biomaterial separation edes], membrane protein crystallization [63], and micro-scale 
gas chromatography [92]. Nanofluidic dynamic arrays have also been devised for high-throughput single nu¬ 
cleotide polymorphism genotyping m- Nanofluidic devices have also been engineered for electronic circuits 
m, local charge inversion [46|, and photonic crystal circuits [35] . Microchannels and micropores have been 
utilized for cell manipulation, cell separation, and cell patterning HUES]. Efforts are given to accomplish 
all steps, including separation, detection and characterization, on a single microchip m- Despite of rapid 
development in nanotechnology, the design and fabrication of nanofluidic systems are essentially empirical at 
present [81]. Since nanofluidic device prototyping and fabrication are technically challenging and financially 
expensive, it is desirable to further advance the field by mathematical/theoretical modeling and simulation. 

The modeling and simulation of nanofluidic systems are of enormous importance and have been a growing 
field of research in the past decade. When the width of a channel is less than 5nm, the transport analysis 
requires the discreteness of substances and, in particular, molecular dynamics (MD) is a useful tool in this 
respect [26]. Typically, the MD determines the motion of each atom in a system using the Newton’s classical 
equations of motion [68 . A simplified model is Brownian dynamics (BD), in which the solvent water molecules 
are treated implicitly, so this method costs less computationally than the MD and is able to reach the time 
scale of physical transport EH EH EH- The BD describes the motion of each ion under frictional, stochastic 
and systematic forces by means of Langevin equation EH EH]. Further reduction in the computational cost 
leads to the Poisson-Nernst-Planck (PNP) theory, which is the most renowned model for charge transport 

The PNP model describes the solvent water molecule as a dielectric 
continuum, treats ion species by continuum density distributions and, in principle, retains the discrete atomic 
detail and/or charge distribution of the channel or pore [6, [33j [59] 104 , 11105 1. The performance of the PB 
model and the PNP model for the streaming current in silica nanofluidic channels was compared [15] . The 
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Brownian dynamics of ions in the nanopore channel was combined with the continuum PNP model for regions 
away from the nanopore channel [2]. The reader is referred to the literature [261 [33] El G31 CG3 1 HQ5] for a 
comprehensive discussion of the PNP theory. A further simplified model is the Lippmann-Young equation, 
which is able to predict the liquid-solid interface contact angle and interface morphology under an external 
electric field [81] . 

Most microfluidic systems involve fluid flow. If the fluid flow through a microfluidic pore or channel is 
also a concern in the theoretical modeling, coupled PNP and the Navier-Stokes (NS) equations can be utilized 
[T6 l [T7 1 123 | |25 | |53 | [88 ] |89 | |9l | [94 ] [96 ] HQ61HQ8] . These models are able to provide a more detailed description 
of the fluid flow away from the microscale pore or channel, i.e., beyond the Debye screening length. 

Recently, a variety of differential geometry based multiscale models were introduced for charge transport 
[941496] . The differential geometry theory of surface provides a natural means to separate the microscopic 
domain of biomolecules from the macroscopic domain of solvent so that appropriate physical laws are applied 
to appropriate domains. Our variational formulation is able to efficiently bridge macro-micro scales and 
synergically couple macro-micro domains [23- One class of our multiscale models is the combination of 
Laplace-Beltrami equation and Poisson-Kohn-Sham equations for proton transport ii. Another class of 
our multiscale models utilizes Laplace-Beltrami equation and generalized PNP equations for the dynamics and 
transport of ion channels and transmembrane transporters [SUM]. The other class of our multiscale models 
alternate the MD and continuum elasticity (CE) descriptions of the solute molecule, as well as continuum fluid 
mechanics formulation of the solvent mmm- We have proposed the theory of continuum elasticity with 
atomic rigidity (CEWAR) [100] to treat the shear modulus as a continuous function of atomic rigidity so that 
the dynamics complexity of a macromolecular system is separated from its static complexity. As a consequence, 
the time-consuming dynamics is approximated by using the continuum elasticity theory, while the less time- 
consuming static analysis is carried out with an atomic description. Efficient geometric modeling strategies 
associated with differential geometry based multiscale models have been developed in both Lagrangian Eulerian 
[36], 37: and Eulerian representations [99] . 

Nevertheless, in nanofluidic modeling, computation and analysis, there are many standing theoretical 
and technical problems. For example, nanofluidic processes may induce structural modifications and even 
chemical reactions P22 IBS], which are not described in the present nanofluidic simulations. Additionally, 
although the PNP model can incorporate atomic charge details in its pore or channel description, which is 
vital to channel gating and fluid behavior, atomic charge details beyond the coarse description of surface 
charges are usually neglected in most nanofluidic simulations. Moreover, as discussed earlier, Stern layer and 
ion steric effect are significant for the EDL, and are not appropriately described in the conventional PNP 
model. Furthermore, nanofluidic simulations have been hardly performed in 3D realistic settings with physical 
parameters. Consequently, results can only be used for qualitative (i.e., phenomenological) comparison and 
not for quantitative prediction. Finally, the material interface induced jump conditions in the Poisson equation 
are seldom enforced in nanofluidic simulations with realistic geometries. Therefore, it is imperative to address 
these issues in the current nanofluidic modeling and simulation. 

The objective of the present work is to model and analyze realistic nanofluidic channels with atomic charge 
details and introduce second-order convergent numerical methods for nanofluidic problems. We present a 
new variational derivation of the governing PNP type of models without utilizing the differential geometry 
formalism of solvent-solute interfaces. As such a domain characteristic function is introduced to represent the 
given solid-fluid interface. Additionally, we investigate the impact of atomic charge distribution to the fluid 
behavior of a few 3D nanoscale channels. We demonstrate that atomic charges give rise specific and efficient 
control of nanochannel flows. Moreover, we develop a second-order convergent numerical method for solving 
the PNP equations with complex nanochannel geometry and singular charges. Furthermore, the change of the 
distribution in atomic charge distribution is orchestrated with the variation of applied external voltage and 
bulk ion concentration to understand nanofluidic currents. Therefore, we are able to elucidate quantitatively 
the transport phenomena of three types of nano-scaled channels, including a negatively charged channel, 
a bipolar channel and a double-well channel. These flow phenomena are analyzed in terms of electrostatic 
potential profiles, ion concentration distributions and current-volt age characteristics. To ensure computational 
accuracy and efficiency for nanofluidic systems, we construct a second order convergent method to solve 
Poisson-Nernst-Planck equations with dielectric interface and singular charge sources in 3D realistic settings. 

The rest of this paper is organized as follows. Section [2] is devoted to a new variational derivation of PNP 
type of models using a domain characteristic function for nanofluidic simulations. In Section [3] we develop a 
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Dirichlet to Neumann mapping for dealing with charge singularities and the matched and boundary interface 
(MIB) method for material interfaces. These methods are employed to compute the PNP equations with 3D 
irregular channel geometries and singular charges. Section [4] is devoted to validate the present PNP calculation 
with synthetic nanoscale channels. We first test a cylindrical nanochannel with one charged atom at the middle 
of the channel and then examine the channel with eight atomic charges that are placed around the channel. 
Since PNP equations admit no analytical solution in general, we design analytical solutions for a modified PNP 
system which has the same mathematical characteristic as the PNP system. In Section [5j we investigate the 
atomic scale control and regulation of cylindrical nanofluidic systems. Three nanofluidic channels, a negatively 
charged channel, a bipolar channel and a double-well channel, are studied in terms of electrostatic potential 
profile, ion concentration distribution and current. Finally, this paper ends with concluding remarks. 


2. Theoretical Models 


Unlike the charge and material transport in biomolecular systems, the charge and material transport in 
nanofluidic systems induces a negligible reconstruction of the solid-fluid interface compared to the system scale. 
Therefore, instead of using our earlier differential geometric based multiscale models [94H96] which allow the 
modification of the solvent-solute interface, we adopt a fixed solid-fluid interface in the present work. To this 
end, we introduce a domain characteristic function in our variation formulation. 

Let us consider a total computational domain 9 C M 3 . We denote 9 m and 9 S respectively the microscopic 
channel domain and the solution domain. Interface T separates 9 m and 9 S so that 9 m |J T |J 9 S = 9. We 
introduce a characteristic function y(r) : M 3 —>► M 3 such that 9 m = y9 and 9 S = (1 — %)9. Obviously, x and 
(1— x) are the indicators for the channel domain and the solution domain, respectively. Unlike the hypersurface 
function in our earlier differential geometry based multiscale models, the interface is predetermined in the 
present model. In the solution domain 9 S , we seek a continuum description of solvent and ions. In the channel 
or pore domain 9 m , we consider a discrete atomistic description. A basic setting of our model can be found 
in Fig. [1} 




(b) 


Figure 1: Illustration of computational geometry, (a) A 3D view of a schematic cylindrical nanochannel whose ends 
are connected to two reservoirs of KC1 solution; (b) A 2D cross-section view of the cylindrical channel whose diameter 
is 10A and length is 49A in the xz- plane. Here, and respectively, represent the applied potential at the left 
end and the right end, and Co represents the bulk ion concentration of both K + and Cl - . 


2.1. Generalized Poisson-Nernst-Planck theory 

Although the PNP theory is quite standard [6j HU H7| |23l [25j |53l [85j [HH1 EH Ell EU 110611108s . it does 
not include non-electrostatic interactions. Here we present a generalized PNP theory by incorporating non¬ 
electrostatic interactions between the solution and the nanoscale channel pore, and between solvent molecules, 
i.e., waters and ions. We utilize a variational formulation to derive generalized PNP equations. 

2.1.1. Energy functional 

Electrostatic energy functional. Electrostatic interactions are ubiquitous at nanoscale and are the dominate 
effects for nanofluidic behaviors. The electrostatic interactions are typically modeled by a number of theoretical 
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approaches, such as the Poisson-Boltzmann (PB) theory [2£l US ESB EZ] , the polarizable continuum theory 
mm and the generalized Born approximation m Eoi- Among these methods, the PB theory is the most 
popular and has a sound origin, i.e., the Maxwell’s equations 0 EE] [70]. A variation formulation of the 
Poisson-Boltzmann theory was originally introduced by Sharp and Honig [76] in 1990 and was extended to an 
electrostatic force derivation [43| and a multiscale formalism [94, 95]. In the present work, we consider the 
following electrostatic energy functional 


C^electr — 



^|V4>| 2 + 4> p m J + (1 — x) 


-||V$| 2 


N c 


$ ^ C o q,;y 


dr, 


(i) 


where 4> is the electrostatic potential, e s and e m are the dielectric constants of the solvent and solute, respec¬ 
tively, and pm represents the fixed charge density of the solute. Specifically, one has p m = Ylk=i Qk$( r — i*/-), 
with Qk denoting the partial charge of the kth atom in the solute and Nf the total number of fixed charges. 
Here C a and g a , respectively, denote the concentration and the charge valence of the <ath solvent species, 
which is zero for an uncharged solvent component. Moreover, N c represents the number of mobile ion species 
through the solution domain. Note that the domain characteristic function x in Eq. 0 is different from the 
hypersurface function S used in our earlier work ]94h96] . 


Non-electrostatic interactions. Non-electrostatic interactions refer to van der Waals interactions, dispersion 
interactions, ion-water dipolar interactions, ion-water cluster formation or dissociation, steric effects, et cetera. 
Some of these interactions are studied in terms of size effects in the past mmmmm m szhbhei. Size 
effects in solvation analysis were accounted with the WCA potential for the solvation p~84f2Q] . Pair particle 
interactions in the Boltzmann kinetic theory and impact to transport equations were formulated by Snider 
et al. in 1996 [79l 180]. To account for solution-channel interactions, as well as ion-water interactions, the 
non-electrostatic interaction energy functional takes the form 

^non- electr L udt =J (1 - x) U dr (2) 


where U is for solvent-channel and ion-ion non-electrostatic interactions. Let assume that the aqueous envi¬ 
ronment has multiple species labelled by a and their interactions with each solute atom near the interface can 
be given by 


Oi 

(3) 

Yi U ak(r) + Y U ap( r )’ 

(4) 


k p 


where C a ( r) is the density of ath solvent component, which may be either charged or uncharged, and U a k is 
an interaction potential between the kth. atom of the channel molecule and the arth component of the solvent. 
For water that is free of salt, C a ( r) is the density of the water molecules. Here U a p( r) is a potential for 
solvent-solvent non-electrostatic interactions, including possible ion-water interactions. 

The solvent-solute interactions in solvation analysis have been represented by the Lennard-Jones potential 
[T8H20] . The Weeks-Chandler-Andersen (WCA) decomposition of the Lennard-Jones potential [93] was utilized 
to split the Lennard-Jones potential into attractive and repulsive parts 


t££’ WCA (r) 


KT WCA (*) 


V LJ 

v Oik > 

Vak +£ak, 

0, 


|r - r fc | < <7fc + cf a , 

|r - r fe | > o-fc +cr a , 

0 < |r — r fe | < <j k + <7a, 
|r - r fe | > cr k + a a , 


(5) 

(6) 


where e a k is the well-depth parameter, and a a are respectively the radii of the kth solute atom and the 
<ath solvent component, r denotes a point on the physical space and represents the location of the kth atom 
in the channel. 

The solvent-solvent interaction term U a p( r) in the total interaction potential U a ( r) does not affect the 
derivation and the form of other expressions. More detailed description of U a p( r) for ion channel transport 
can be found in our earlier work [15] [96:. 
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Chemical potential related free energy. Chemical potential related free energy is essential for the description 
of mobile charges in the nanofluidic system 

Gchem = £ 51 { M C » + k sTC a In (PGj - k B T (C a - C a0 ) J dr, (7) 

where is a reference chemical potential of the <ath species at which the associated ion concentration is C a o 
given = U a = /a a o = 0 and /a a o is a relative reference chemical potential which is the difference between 
the equilibrium concentrations of different solvent species. Here ks is the Boltzmann constant and T is the 
temperature. The term ksTC a In (^§^j is the entropy of mixing, and —ksT (C a — C a o) is a relative osmotic 
term [66] . 

It is standard to determine the chemical potential of species a by the variation with respect to C a [96 

^G c hem chem 0 , 7 m 1 ^ol / c \ 

—^—=^/i a -/iaO + ^Tln—. ( 8 ) 

Note that at equilibrium, /i^ hem 7 ^ 0 and C a 7 ^ C a 0 because of possible external electrical potentials, solvent- 
solute interactions, and charged species. Even if the external electrical potential is absent and system is at 
equilibrium, the charged solute may induce the concentration response of ionic species in the solvent so that 

C a ± C a 0- 


Total energy functional. The total free energy functional for the nanofluidic system consists of the electrostatic 
interactions, non-electrostatic interactions and chemical potential related energy term 


G3x,<MO>}] = J lx 


+(1 -x)U 


-^|V $| 2 + $p m j+(l-x) 


c„ 


N c 


y |V$| 2 + $ 55 C aQa 


(9) 


+(1 - X) 5Z « - M«o) c a + k B TC a In - k B T (Ca - C a0 ) + A a C, 


• dr, 


where the first row is the electrostatic free energy functional, the second row is the free energy functional of 
non-electrostatic interactions, and the third row is chemical potential related energy functional. Note that the 
electrostatic free energy functional is the same as the polar solvation free energy functional QUES US]. Here 
A a is a Lagrange multiplier, which is included to ensure appropriate physical properties at equilibrium [38] . 

Although it appears that the characteristic function y is quite similar to the hypersurface function in our 
earlier work [94H96] . x is j us t an indicator for a given molecular domain with a fixed interface T. In 
contrast, the hypersurface function in our earlier work not only as a characteristic function for the molecular 
domain, but also plays the role of a moving interface whose dynamics is driven by mechanical and electrostatic 
forces, i.e., the Laplace-Beltrami equation. However, the fixed interface T in the present work is a reasonable 
approximation for nanochannels. 

An important feature of the present total free energy function formulation © is that the free energy 
functional U of the non-electrostatic interactions is employed for nanofluidics. Therefore, the present theory 
is able to deal with a variety of non-electrostatic interactions. Consequently, the solvent microstructure near 
the channel can be predicted correctly. 


2.1.2. Governing equations 

The total free energy functional (|9| is a function of electrostatic potential $ and ion concentration C a - 
Unlike our earlier formulations [94H96] T the solvent-channel interface T is given in the present work. Like our 
earlier work, governing equations for the system are derived by using the variational principle. 

The Poisson equation. The variation of the total free energy functional with respect to the electrostatic 
potential <f> results in the classical Poisson equation 

N c 

~ V- (e(x)V$) =XPm + a-x)T, Caqa ’ r e ( 10 ) 

OL 
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where e(x) = (1 — x) e s -\~x e m is the dielectric profile, which is e m in the molecular domain and e s in the solvent 
domain. Due to the characteristic function x, the Poisson equation (10) can be split into two equations 


- v • (e m V$) = pm, rG n m (11) 

N c 

—V • (e s V$) = Y^C a q a , r£(l s . (12) 

a 

However, electrostatic potential <h(r) is defined on the whole computational domain (for all r G D), Therefore, 
at the solution-channel interface T, the following interface jump conditions are to be implemented to ensure 
the well-posedness of the generalized Poisson equation Ha uni Ena 

Mr)] = 0, r e T (13) 

[e(r)V$(r)] • n = 0, r G T (14) 


where [•] denotes the difference of the quantity cross the interface T. n is the unit norm of T and 

<«=< i5 > 


In Eq. (10), the densities of ions C a are to be determined by the variational principle as follows. The boundary 


conditions of Eq. (10) depend on experimental settings. Typically, mixed boundary conditions are used. 


Generalized Nernst-Planck equation. It is also standard to determine the relative generalized potential /i| en 
by the variation with respect to the ion density C a 


C 

= l 1 a ~ Ah*o + A^Tln -P- + q a <f> + U a + A a = // 

^c*0 


chem 


+ + U a + A 0 


(16) 


At equilibrium, we require /i| en , rather than /i^ hem , to vanish 


^ oc — 

c a = C a Qe 


Qa^ + Ug-^a 0 

k B T 


(17) 


Therefore, the relative generalized potential /i| en is simplified as 

/x® en = k B T\n + q a <t> + U a - HaO- 


(18) 


We derived a similar quantity from a slightly different perspective in our earlier work Note that the 

relative generalized potential consists of contributions from the entropy of mixing, electrostatic potential, 
solvent-solute interaction and the relative reference chemical potential. In practice, the nanofluidic system is 
out of equilibrium due to applied external field and/or inhomogeneous concentration across the nanochannel. 
Therefore, /a^ n does not vanish. In general, a major mechanism for establishing the equilibrium is diffusion 
processes driven by gradients of density, velocity, temperature and electrostatic potential [79l [80] . By Fick’s 

g en 

first law, the ion flux of diffusion can be given as J a = —D a C a \/^p^ with D a being the diffusion coefficient 
of species a. We therefore have the diffusion equation for the mass conservation of species a at the absence 
of steam velocity 


dC a 

dt 


V • J a . 


In the explicit form, the generalized Nernst-Planck equation is 


dCa 

dt 


= V- 


D a 


k B T 


V {q a $ + U a ) 


a =»■!,•■ 


,#c, 


(19) 


where q a + U a can be regarded as the potential of the mean field. At the absence of non-electrostatic 

Lees to the standard Nernst-Planck equation. At the steady state, one has 

D a (\7C a + kh;V(q a <S> + U a ) )|=0, a = l,---,N c . (20) 




k B T 
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Note that Eq. (19) does not involve the characteristic function x because it has already been restricted 
to solution domain Q s . In contrast, generalized Poisson equation (10) is defined on the whole domain (D). 
The nature boundary condition is assumed in our derivation. However, due to the experimental setup, mixed 
boundary conditions are typically employed in our simulations. 


3. Mathematical Algorithms 

The geometric setting of the ionic diffusive nanofluidic system employed in the present investigation is 
described. In this work, we develop a second-order PNP solver for 3D ionic diffusive nanofluidic channels 
with irregular geometry and material interface. To this end, we appropriately modify the computational 
algorithms developed in our earlier work nna for nanofluidic systems. To emphasize the primary effects of 
atomic charges in ion diffusive nanofluidic channel design and the use of interface techniques in 3D nanofluidic 
systems, we neglect the non-electrostatic interactions, i.e., assuming U = 0 , in the rest of this work. However, 
since non-electrostatic interactions are important for nanofluidic systems mm, the situation that U ^ 0 will 
be investigated in our future work. To avoid confusion, our generalized PNP model works for all kinds of 
ion species. However, our model is designed based on realistic transmembrane channels. Since potassium 
and chloride are most important ion species in cellular charge transport, we use the potassium chloride (KC1) 
system as an example in our model. 

3.1. A schematic diagram of ionic diffusive nanofluidic channels 

In the present work, we construct a cuboid nanofluidic system whose dimensions are 16 x 16 x 56A 3 . It 
contains a cylindrical nanochannel which is placed at the center of the system as illustrated in Fig. 0a). The 
radius of the channel pore is 5A and the length of the channel is 49A as depicted in Fig. 0 b) . Each of 
the channel ends is connected to a reservoir of potassium chloride (KC1) solution. In our simulation, the 
computational domain Q is the nanofluidic system and it is mainly divided into ion inclusion region Q s and 
ion exclusion region . The ion inclusion region is the region inside the channel and the two reservoirs where 
ions may penetrate and travel through. The ion exclusion region is the rest where there is no mobile ion, 
but has fixed charged particles. In contrast to our differential geometry based multiscale models [941496] . the 
interface between two regions $d s and is denoted by T, which is predetermined by the channel structure 
and does not change during our simulation. 

A number of properly charged atoms about 1.8A apart are positioned around the channel so that the 
channel flow can be regulated by electrical charges. In reality, these charged atoms can be realized by appro¬ 
priate dopants. The ^-coordinate of the atoms along the channel length is determined first and then at each 
cross section perpendicular to the channel axis, the atoms are aligned along a concentric circle whose size is 
sufficiently bigger than that of the channel pore. The locations of the atoms are equally spaced according to 
the circumference of the circle. Figure [2] shows an example of placing four negative charges around the channel 
at z = 0 A. In the cross section on the xy- plane, we divide a concentric circle with radius 6.5A into four parts 
and then locate each anion as described in Fig. |2|a). Managing the number, magnitudes and signs of charges 
enables one to generate various types of atomic charge distributions for the cylindrical nanofluidic channel. 


3.2. Iterations of Poisson and Nernst-Planck equations 

The MIB method is utilized to solve the interface problems Eqs. (10) and (19). Since these equations 
are coupled, an iterative procedure is required to obtain convergent results. Here we outline this solution 
procedure. To ensure that the iteration is convergent, 4> and C a are updated by a successive over relaxation 
(SOR)-like iterations 

j $ new = (1 - w p ) $ new + w p $ old 
{ C" ew = (1 - W C ) C" ew + w c C° ld , 


where w p and w c are appropriately selected between 0 and 1. This iterative procedure is efficient for ion 
channel problems [18, 104 . The relaxation parameters w p and w c should be sufficiently close to 1 to ensure 
the convergence. They have little influence on computational results as long as the iteration is convergent m- 
Although the change of applied voltage or bulk ion concentration may requires a number of iteration steps, 
the convergence is still maintained m- In this computation, the values for both relaxation parameters are 
fixed at w p = w c = 0.9. If the iteration is divergent, we adjust these values to w p = w c = 0.8. 
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(a) (b) 

Figure 2: Illustration of charge distribution, (a) Four negatively charged atoms are equally placed around the 
cylindrical nanochannel; (b) A 3D schematic diagram consisting of the cylindrical channel with four charged atoms. 


After the electrostatic potential and the ionic concentration C a are iteratively solved, the electric current 
is computed at each cross section inside the channel along the channel axis m- For each ionic species <a, its 
current is calculated by 



d_Cg 

dz 


q a r 
k B T a dz 


Ax Ay, 


( 21 ) 


where S is the cross section in the xy- plane. The total current is the sum of two ionic currents. Actually, 
there is no significant change along the location of the cross section and hence the current through the center 
of the channel axis is usually chosen to elucidate the current-voltage (I-V) relation. 


4. Convergent Validation 

In this section, we construct analytically solvable systems to validate the proposed numerical methods. 
The analytic solution of the PNP equations is unknown for realistic geometries. However, it is a standard 
procedure to design analytical solution for slightly modified PNP equations which share the same mathematical 
characteristic with the original PNP equations [ 104 . Consequently, the numerical convergence of designed 
solution algorithms can be validated. 

We first present the analytical solution to a set of PNP-like equations. Additionally, we consider two simple 
examples, one with a single atomic charge, and the other with eight atomic charges, to verify the second order 
convergence of our numerical methods. Finally, both a negatively charged nanofluidic channel and a bipolar 
nanofluidic channel are utilized to further validate the proposed numerical methods. 

We set e m = 1 and e s = 80 in all the numerical tests in this section. Therefore, there is a sharp discontinuity 
in the dielectric coefficients across the solvent-solute interface. 


4.1. Analytical solution system 

We consider a set of Nf charged atoms at with fixed charge Q k , where k = 1,2, •• • ,7Vy in Q m . The 
geometry of the analytically solvable system can be arbitrary in principle. However, one can refer to the 
cylindrical geometry described in Section [3T] with the cross section of the cylinder illustrated in Fig. [3J We 
set the electrostatic potential to 


<£> (r) = < 


N f 


Qk 


K 3e 




cos x cos y cos z + ^_ 

*=i (x - x k f + (y - y k f + 0 - z k f 

0 . 4 ^ 

cos x cos y cos z, r E il s , 


( 22 ) 
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Figure 3: A cross section of the cylindrical geometry of the nanofluidic channel. The ion inclusion region is 
composed of two reservoirs and the inside of the cylindrical channel. Red circles indicate the charged atoms around 
the channel to generate channel surface charge effects. 


where r = (x,y,z) and (xk,yk, z k) is the atomic central position of the kth charge near the channel surface. 
The concentrations for two mobile ion species are set to 


C\ (r) = 0.2 cos x cos y cos 2: + 0.1 
62 ( 1 *) = 0.1 cos x cos y cos 2: + 0.1 


only for r E £2 S , 


(23) 


and Ci(r) = 62 ( 1 *) = 0 for all r E because ions are able to move only in the solution confined in the 
region Q s . Indeed, this set of solutions satisfies the following PNP-like equations 

N f 

-V • (e (r) V3> (r)) = 4ir ^ Q k 6 (r - r k ) + [C\ (r) - C 2 (r)] + R (r) in Q, 

< k =1 (24) 

V • [VCi (r) + Cl (r) V$ (r)] = R x (r) in 

k v • [VC 2 (r) - C 2 (r) (r)] = R 2 (r) in tl s , 


where we have set qi = 1, #2 = — 1, and Di( r) = D 2 (r) = 1. Moreover, for every r E 


' R( r) = —3 cos x cos y cos z 

< Ri (r) = —0.6 cos x cos y cos z + V • [(0.2 cos x cos y cos z + 0.1) V (cos x cos y cos z)\ 

^ R 2 (r) = —0.3 cos x cos y cos z — V • [(0.1 cos x cos y cos z + 0.1) V (cos x cos y cos z)\ . 

However, R( r) = Ri(r) = R 2 (r) = 0 for all r E Q m . Additionally, the jump conditions at the interface T 
can be specifically given as the follows: 

[<h (r)] = <f>* (r) + fl — cos x cos y cos 2: 

J V s / 

[e (r) ^> n (r)] = e m V [cos x cos y cos z + <h* (r)] • n + e s V cos x cos y cos z'j • n, 

where n is the outward unit normal vector. 
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In order to investigate the convergence order, we apply two error measurements 


L 


oo 


max I F, 

i,j,k 


num 

i,j,k 


7T , exact 


and L 2 



where FfYk and respectively, represent the numerical and exact values of a function F at (xi,yj,Zk) 

and N is the total number of computational nodes. 


4-2. A cylindrical nanochannel with a single atomic charge 



Figure 4: Illustration of the geometries of two test cases, (a) A cylindrical nanochannel with a single negative charge; 
(b) A cylindrical nanochannel with eight negative charges. 


We first test the cylindrical channel with a single atomic charge which is placed at (6.5, 0, 0) and whose 
charge is —0.08e c as shown in Fig. Qa). The set of analytical solutions of the PNP equations introduced in 
Section 4.1 is used to compare with numerical results by solving Eq. (??). Table [l] demonstrates numerical 
errors and convergence orders for different number of computational nodes. The fixed charge of the atom 
influences on the accuracy of the electrostatic potential 4> and the negative charge of the atom enhances the 
errors and reduces the orders of the negative ion concentration C 2 - It is interesting to note that the simulation 
attains a good second order convergence. 


4-3. A cylindrical nanochannel with eight atomic charges 

Next, we explore this cylindrical channel with eight atomic charges outside the nanochannel. Consider two 
cross-sections perpendicular to the channel length at 2 = — llA and z = llA as illustrated in Fig. (4^b). In 
order to put negative ions same distance away from the cylinder surface and same angle difference between 
two atoms, we employ the polar coordinate system. The distance between each atom and the origin is always 
6.5A and the angle from the positive x-axis is increased by a right angle. Therefore, the coordinates of these 
atoms are as follows: 

^6.5cos ^(i — 1)^ ,6.5sin ^(i — 1)^ , —11^ and ^6.5cos — 1)^ ,6.5sin — 1)^ , 11^ 

for i m 1, 2, 3, and 4. At each cross section, we obtain four point charges that are equally spaced. As shown 
in Table [2j the errors and orders in solving the PNP equations with this atomic charge setting generates 
little difference from those with a single atomic charge. This validation test also indicates the second order 
convergence of our methods. 
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Table 1: Numerical errors and orders for the cylindrical channel with a single atomic charge. 



Mesh size 

Too 

L 2 

Error 

Order 

Error 

Order 

$ 

II 

o 

1.3742E-2 

- 

3.1647E-3 

- 

h = 0.32 

8.7238E-3 

2.0362 

2.0215E-3 

2.0088 

h = 0.2 

3.5614E-3 

1.9062 

8.7926E-4 

1.7712 

h = 0.16 

2.2282E-3 

2.1016 

5.3548E-4 

2.2224 

Cl 

II 

o 

3.8661E-3 

- 

1.0406E-3 

- 

h = 0.32 

2.5648E-3 

1.8390 

6.4972E-4 

2.1106 

h = 0.2 

8.4896E-4 

2.3524 

2.5014E-4 

2.0309 

h = 0.16 

5.4286E-4 

2.0039 

1.5920E-4 

2.0250 

c 2 

II 

o 

2.0352E-3 

- 

5.3493E-4 

- 

h = 0.32 

1.3130E-3 

1.9642 

3.2772E-4 

2.1958 

h = 0.2 

4.7824E-4 

2.1488 

1.3661E-4 

1.8617 

h = 0.16 

2.7990E-4 

2.4327 

8,0089E-5 

2.3930 


Table 2: Numerical errors and orders for the cylindrical channel with eight atomic charges. 



Mesh size 

Tqo 

L 2 

Error 

Order 

Error 

Order 

<f> 

II 

o 

1.3745E-2 

- 

3.1649E-3 

- 

h = 0.32 

8.7249E-3 

2.0369 

2.0214E-3 

2.0092 

h = 0.2 

3.5571E-3 

1.9090 

8.7894E-4 

1.7720 

h = 0.16 

2.2300E-3 

2.0925 

5.3564E-4 

2.2194 

Cl 

II 

o 

3.8661E-3 

- 

1.0406E-3 

- 

h = 0.32 

2.5648E-3 

1.8390 

6.4972E-4 

2.1106 

h = 0.2 

8.4892E-4 

2.3525 

2.5015E-4 

2.0308 

h = 0.16 

5.4286E-4 

2.0037 

1.5920E-4 

2.0252 

c 2 

II 

o 

2.0353E-3 

- 

5.3493E-4 

- 

h = 0.32 

1.3130E-3 

1.9643 

3.2772E-4 

2.1958 

ft = 0.2 

4.7824E-4 

2.1488 

1.3659E-4 

1.8620 

h = 0.16 

2.7790E-4 

2.4326 

8.0091E-5 

2.3923 


4-4- A negatively charged nanofluidic channel 

Now, we perform another numerical test to verify the convergence and accuracy of the proposed PNP cal¬ 
culation on nanofluidic channels. A negatively charged nanochannel, or a unipolar nanochannel is constructed. 
The channel length on the z-axis is divided into 27 subdivisions. At each circular cross section, eight charged 
atoms which are 1.5A apart from the cylinder surface are equally spaced. Each atom has a charge of — 0.08e c 
and is about 1 . 8 A apart from other charges. 

First of all, we examine colored surface plot and contour plots of electrostatic potential distributions 
along the negatively charged channel by computing Eq. (??). The computational results are demonstrated in 
Fig.|5ja) and Fig. [ 6 j Figure[5]is useful to understand the atomic charges of the nano-scaled channel. Moreover, 
it is interesting to notice that our proposed PNP solver works very well with nanofluidic channels from the 
fact that the boundary line of the channel are obvious to identify in contour plots. 

Figure |5ja) shows the electrostatic potential profiles over the surface of the negatively charged channel. 
Most of the parts of the channel surface have blue colors. Since the blue colors indicate the negative electrostatic 
potential values, it is obvious that the channel surface possesses negative charge. Additionally, three contour 
plots are described in Fig. [ 6 ] at ^ = —10A, z = OA and z = 10A. In every picture, the right outside of the 
channel is dark blue, which also implies the channel surface is negatively charged. 

Then we use the same analytical solutions (??) and (??) to find the numerical errors and orders. The results 
are given in Table [3j Through this test, we observe that our proposed PNP numerical schemes achieve the 
second order accuracy in computing the potential and ion concentrations for the negatively charged channel. 

4-5. A bipolar nanofluidic channel 

We also consider a bipolar nanofluidic channel which functions as a nanofluidic diode. It is a nano-sized 
channel whose atomic coordinates are the same as those of the aforementioned negatively charged nanochannel. 
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(a) A negatively charged channel (b) A bipolar channel 

Figure 5: Illustrations of computational results of electrostatic potential profiles on channel surfaces. Here, the blue 
colors represent negative values, while the red colors represent positive values. The negatively charged channel surface 
is mostly blue, but the color of the bipolar channel surface is changed from red to blue. 
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c 

< 


if) 

O) 
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< 


x (Angstrom) 

(b) Cross section: 2 = OA 

Figure 6: Three contour plots of electrostatic potential profiles of the negatively charged channel at z = — 10A, 
z — OA and z — 10A on the xy-plane when the mesh size is h = 0.2A. 


x (Angstrom) 

(a) Cross section: 2 = —10A 


x (Angstrom) 

(c) Cross section: z = 10A 
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c 
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However, its charges are altered from positive to negative or vice versa at the middle of the channel axis mm- 
In our bipolar channel, the first half cylinder is affected by atoms with charge of 0.08e c and the atomic charges 
on the other half are —0.08e c . 

The computational results of the electrostatic potential through the bipolar channel is shown in Fig. [5jb) 
and Fig. [7J From Fig. |5jb), we are able to see that the atomic charges are changed from positive to negative. 
Such properties of the bipolar channel are also clearly manifested in the cross-sectional results in Fig. [7| When 
we compare the contour plots in Fig. [7^a) and Fig. [7^c), it is a little bit difficult to distinguish the colors of 
the outside of the bipolar channel. However, the channel inside clearly shows different colors, which indicates 
that the change of atomic charges influences the ion transport through the channel. 

The validation of our numerical methods for the bipolar channel is given in Table [4] Again, we see a good 
second-order convergence for this test problem. In the next section, we apply our PNP simulator to study 
three nanochannels. 


14 





















Table 3: Numerical errors and orders for the negatively charged channel. 



Mesh size 

Loo 

L 2 

Error 

Order 

Error 

Order 

$ 

II 

o 

1.3803E-2 

- 

3.1672E-3 

- 

h = 0.32 

8.7573E-3 

2.0389 

2.0218E-3 

2.0116 

h = 0.2 

3.5184E-3 

1.9401 

8.7433E-4 

1.7836 

h = 0.16 

2.3190E-3 

1.8682 

5.3538E-4 

2.1980 

Cl 

II 

o 

3.8725E-3 

- 

1.0406E-3 

- 

h = 0.32 

2.5664E-3 

1.8437 

6.4976E-4 

2.1103 

ft = 0.2 

8.6571E-4 

2.3121 

2.5037E-4 

2.0291 

h = 0.16 

5.4302E-4 

2.0902 

1.5936E-4 

2.0244 

c 2 

II 

o 

4^ 

2.0402E-3 

- 

5.3508E-4 

- 

h = 0.32 

1.3126E-3 

1.9765 

3.2767E-4 

2.1978 

h = 0.2 

4.8312E-4 

2.1265 

1.3660E-4 

1.8615 

h = 0.16 

2.7784E-4 

2.4793 

8.0291E-5 

2.3816 



-8 -6 -4-2 0 2 

x (Angstrom) 


• 

8 

6 

• •••• 


? 4 
2 2 



• 

+■* 

V) n 
O) U 

c _ 

• 

• 

• 

<-2 

>- 4 


jm 

• 

-6 

-ft 

• •••• 


? 4 
2 2 

</> n 
O) U 
C 

<-2 
>_4 
-6 


-8 -6 -4 -2 0 2 4 6 8 

x (Angstrom) 


- 8 , 



I 


-0.5 


-1 


I 


-8 -6 -4 -2 0 2 4 6 8 


(a) Cross section: z = — 10A (b) Cross section: z = 0A (c) Cross section: z = 10A 

Figure 7: Three contour plots of electrostatic potential profile through the bipolar channel at 2 = —10A, z = 0A and 
2 = 10A on the xy- plane when when the mesh size is h — 0.2A. Here, the blue colors represent negative values, while 
the red colors represent some small positive values. At the cross section of z = — 10A, the inside of the channel is blue 
(i.e., negatively charged); on the other hand, at the cross section of z — 10A, the inside of the channel is red (i.e., 
slightly positively charged). 


Table 4: Numerical errors and orders for the bipolar channel. 



Mesh size 

Loo 

l 2 

Error 

Order 

Error 

Order 

$ 

II 

o 

4^ 

1.3752E-2 

- 

3.1649E-3 

- 

h = 0.32 

8.7568E-3 

2.0227 

2.0217E-3 

2.0086 

h = 0.2 

3.5649E-3 

1.9121 

8.7942E-4 

1.7711 

h = 0.16 

2.2421E-3 

2.0782 

5.3585E-4 

2.2201 

Ci 

II 

o 

3.8743E-3 

- 

1.0406E-3 

- 

h = 0.32 

2.5663E-3 

1.8458 

6.4975E-4 

2.1105 

h = 0.2 

8.5125E-4 

2.3479 

2.5016E-4 

2.0308 

h = 0.16 

5.4306E-4 

2.0143 

2.5016E-4 

2.0191 

c 2 

II 

o 

2.0384E-3 

- 

5.3493E-4 

- 

h = 0.32 

1.3132E-3 

1.9706 

3.2773E-4 

2.1967 

h = 0.2 

4.8477E-4 

2.1203 

1.3665E-4 

1.8611 

h = 0.16 

2.7790E-4 

2.4935 

8.0280E-5 

2.3837 


5. Simulation Results 

Having validated the numerical convergence of our proposed PNP algorithm, we explore the nature of 
charged nanofluidic channels under various physical conditions including applied voltage, atomic charge distri- 
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(a) (b) 

Figure 8: Effect of applied voltage on a negatively charged nanofluidic channel, (a) Electrostatic potential profiles; 
(b) ion concentration distributions along the channel length (z— axis). Here A<f> is varied from OV (square), 0.2V 
(triangle), 0.4V (asterisk), 0.6V (diamond), 0.8V (circle) to IV (plus sign), where A4> is the difference of the applied 
voltage between two ends of the system. The charge of each atom near the channel surface is Qk = —0.08e c and bulk 
concentration is Co = 0.01M for both ions K + and Cl - . Two dashed vertical lines indicate the ends of the cylindrical 
nanochannel. As the applied voltage difference A<f> gets larger, the potential at the right part of the inner channel is 
increased. Consequently, more K + ions (dashed line) are accumulated at the left part of the inner channel. However, 
there is little change in the concentration of Cl - ion (solid line). 


bution and bulk ion concentration. We investigate ion concentration distributions and electrostatic potential 
profiles along the channel direction (^-direction) and their values are averaged over the xy -cross section at 
each z. We also illustrate current-volt age (I-V) curves in which the current at the center of the channel pore is 
used. Particularly, ion concentration distribution describes the movements of two different ion species through 
a channel in detail and current-voltage characteristic clearly shows electrical features of a nanochannel. We 
examine three kinds of nano-scaled channels, namely a negatively charged channel, a bipolar channel and a 
double-well channel. We have used e m = 2 and e s = 80 for all computations in this section. 

5.1. A negatively charged nanofluidic channel 

In order to clarify the role of atomic charges in nanofluidic systems, we first consider a negatively charged 

The atomic charges of the channel are specified in Section |4.4| Here, all 
of the ions around the channel have the negative sign. We vary the voltage at the end of the right reservoir 
and keep bulk ion concentrations of both reservoirs unchanged. Therefore, applied voltage is the driving force 
to relocate potassium ions and chloride ions within the nanochannel. The atomic charges determine the ion 
selectivity of the nanochannel and bulk ion concentration affects the current migrated through the channel. 

5.1.1. Effect of applied voltage 

First, we study the impact of applied external voltage to behavior of ions traveling within a negatively 
charged nanochannel. We set the bulk ion concentrations of K + and Cl - as Co = 0.01M and the charge of 
each atom placed around the channel surface as Qk = —0.08e c . The voltage applied at left end of the system, 
4V, is fixed at 0V and the voltage applied at right end of the system, 4V, is increased gradually from 0V to 
IV. The A4> represents the difference between 4V and 4V, where A<f> = 4V — 4V- Figure [8] illustrates the 
electrostatic potential and ion concentrations along the z—axis at the center of the channel pore in response to 
the external voltage difference. In Fig.[8ja), as 4V gets increased, the electrostatic potential at the right part 
of the inner channel becomes dramatically higher. As a result, more cations are accumulated at the left part 
of the inner channel as shown in Fig.^b). In fact, the negative atomic charges of the channel electrostatically 
repel anions and attract cations, which makes the solution within the channel almost unipolar. 


channel described in Section 4.4 
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Figure 9: Ionic current for each ion species versus the applied potential difference A<f>. (a) The transport behavior of 
the nanochannel without atomic charge; It should be noticed that we use different relative reference chemical potential 
/la o for different ion species, Therefore there is a small current difference between potassium and chloride ions, (b) 
The transport behavior of a negatively charged channel when Qk — — 0.08e c . Here, the bulk ion concentrations at 
both reservoirs Co = 0.01M are fixed. When Qk = 0e c , both current-voltage graphs are linear and the positive current 
is roughly double of the negative current; on the other hand, when Qk = —0.08e c , the positive current-voltage graph 
(square) is nonlinear and the negative current-voltage graph (triangle) is almost always zero. Moreover, the positive 
current is much larger than the negative current and the difference gets increased as the voltage increases. 


By comparing the current-voltage (I-V) curve of the negatively charged channel with that of an uncharged 
one, we are able to clarify the effect of atomic charges in a nanometer channel. In these two graphs, the 
current values are obtained using Eq. (??) at the cross section inside the channel which lies in the xy -plane 
when z = 0A. Figure [9|a) gives the relationship between current and voltage of each ion species for the same 
dimensional nanochannel without atomic charge (Qk = 0e c ). Both of the ionic currents are proportional to the 
applied voltage and the K + current is roughly double of the Cl - current. In this case, the nanochannel obeys 
the Ohm’s law and is non-selective. However, the negative atomic charges destroy the linearity of the positive 
current-volt age characteristics and generate a large difference between two ionic currents as shown in Fig.^b). 
This nonlinearity or deviation from the Ohm’s law is closely related to the non-proportionality between the 
potential change within the channel and the applied voltage [28] . Since the negative atomic charges near the 
channel surface hinder the access of Cl - ions, the negative current is almost zero for every applied voltage. 
Thus we can conclude that a nanochannel with unipolar atomic charges generates a charged current mostly 
composing of counterions which can be increased by providing more external voltage. 

5.1.2. Effect of atomic charges near the channel surface 

Next, we examine the effect of atomic charges on ionic flow through the negatively charged channel. Except 
for the magnitude of Qk, we fix the bulk ion concentrations for both ion species at Co = 0.01M and the applied 
voltage difference at Aff> = IV. A stronger negative atomic charge, i.e., a larger value of \Qk\, encourages the 
cation accumulation inside the channel and prevents the anions from entering the channel. As a result, the 
K + current is increased and the Cl - current is decreased to near zero. In fact, the positive current is not 
directly proportional to the magnitude of atomic charge because of a stronger diffusion induced by the larger 
concentration gradient [28 . It is interesting to note that the channel current can be controlled by the atomic 
charge amplitude. When the amplitude reaches a suitable threshold, almost all ions with the same sign of 
charge with the channel atomic charge cannot penetrate through the inlet of the nanochannel. Therefore, the 
proposed nanochannel has a near perfect ion selectivity. 
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Figure 10: Effect of atomic charges on a negatively charged nanochannel. The ionic current in response to the change of 
the charge Qk of the atoms placed around the channel surface is depicted. Six different charges Qk — 0e c , Qk = —0.02e c , 
Qk — —0.04e c , Qk — —0.06e c , Qk — —0.08e c and Qk — —0.1e c are simulated when the bulk concentration is Co = 0.01M 
and the applied voltage difference is A<f> = IV. Here, \Qk \ is the magnitude of the charge of the atoms. As the magnitude 
of the atomic charges is increased, the ionic current of K + is sharply amplified as well. However, the ionic current of 
Cl - is reduced to near zero. 



Figure 11: Effect of bulk ion concentration on a negatively charged nanochannel. The total channel current versus 
the external voltage difference (I-V) is shown at five different bulk concentrations Co = 0.01M (square), Co = 0.05M 
(triangle), Co — 0.1M (diamond), Co — 0.2M (circle) and Co — 0.4M (plus sign) when Qk = — 0.08e c . As the bulk ion 
concentration is increased, the total current gets higher and the I-V characteristics becomes near linear. 
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5.1.3. Effect of bulk ion concentration 

Another important aspect which is necessary to understand the transport within a nanofluidic channel is the 
bulk ion concentration. The electrical double layer produces a unique difference between a nanofluidic channel 
and a microfluidic one. The bulk ion concentration is a crucial factor to determine the thickness of the double 
layer. In fact, when double layers overlap inside a nanochannel, the aqueous solution confined in the channel 
expresses charge opposite to the atomic charges of the nanochannel [SBj. Figure 11 shows the total current as 
a function of the applied voltage difference for five different bulk ion concentrations Co = 0.01M, Co = 0.05M, 
Co = 0.1M, Co = 0.2M and Co = 0.4M when the atomic charge Qk is assumed to be — 0.08e c . As the bulk 
ion concentration is increased, more cations penetrate through the channel and thus the total channel current 
is dramatically increased. A higher bulk ion concentration reduces the double layer and the channel surface 
becomes neutralized by the pulled counterions within the layer [28]. Consequently, the I-V graph becomes near 
linear, i.e., obeying Ohm’s law. It is noted that a charged channel with high bulk ion concentration behaves 
like an uncharged one, due to the decrease in the Debye length. This phenomenon is obviously manifested in 
the positive ion concentration distributions in Fig. [l2j Figure [l2ja) shows the concentration profiles of the K + 
ion along the x-direction at four different bulk ion concentrations Co = 0.01M, Co = 0.05M, Co = 0.1M and 
Co = 0.2M. Figure (I^b) demonstrates the normalized concentration profiles, which are generated by scaling 
each ion concentration with its bulk value C(r)/Co- As illustrated in the figure, the Debye screening effect 
can be observed by the concentration distributions, which are higher at the channel boundary and lower at 
the channel center. Additionally, for the normalized concentration profiles, the lower the bulk concentration, 
the higher the normalized concentration, which indicates the larger Debye-Layer. 



(a) 



x (Angstrom) 

(b) 


Figure 12: The effect of bulk ion concentration on the double layer effect, (a) Ionic concentration distributions 
of the positive ion along the x-axis at four different bulk ion concentrations Co = 0.01M, Co — 0.05M, Co = 0.1M 
and Co = 0.2M with Qk = — 0.08e c and AT = 0.4V being fixed, (b) Normalized ionic concentration distributions of 
(a)indicates the large Debye-Layers for lower bulk ionic concentrations. 


The normalized current is considered to assure that a higher bulk ion concentration weakens the role of 
atomic charges. The normalized current is computed by dividing the current through the negatively charged 
channel by that through the uncharged one at several bulk ion concentrations. Excluding the atomic charges, 
both channels have the same voltage difference and the same bulk ion concentration. Figure [T3^a) represents 
the normalized ionic current for each ion species and Fig. [I3|fo) shows the normalized total channel current. 
In both figures, the normalized values get close to one as the bulk ion concentration is increased, which implies 
that the charged channel with a high bulk ion concentration does not demonstrate much atomic charge effect 
in the transport phenomena. 
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Figure 13: The effect of bulk ion concentration on the normalized channel current for a negatively charged nanochan- 
nel. (a) The normalized ionic current; (b) The normalized channel current with respect to the increase in bulk ion 
concentration. The atomic charge Qk — — 0.08e c and the applied voltage difference A<f> = IV are fixed. The normalized 
current is the quotient of the current of the negatively charged channel and the current of the uncharged channel when 
Qk = 0e c . As the bulk ion concentration gets larger, the normalized value becomes almost one. The negatively charged 
channel with a sufficiently larger bulk ion concentration behaves like a uncharged channel. 


5.2. A bipolar nanofluidic channel 

In this section, we examine a bipolar channel whose size and atomic charge constitution are described in 
Section [45] In this channel, the first half of the channel is positively charged and the second half is negatively 
charged. A bipolar channel can behave like a p-n junction. Therefore, it is interesting to explore the transport 
properties of the bipolar nanochannel. 


5.2.1. Effect of applied voltage 

We first consider three types of voltage bias across the channel length. One case is called no bias in which 
both ends of the system have zero voltage. Another case is named a forward bias, for which the voltage applied 
at the right end of the system is IV. The other case, on the contrary, is referred a reverse bias, for which the 
voltage at the left end of the system is IV. Additionally, we fix the bulk ion concentration of KC1 at 0.1M and 
set the amplitude of atomic charge \Qk\ to 0.08e c . Figure 14 compares the electrostatic potential profiles and 
ion concentration distributions along the ^-direction in each case. 

As shown in Fig. p!4^i-a), when A<f> = OV, the electrostatic potential is high under the positive atomic 
charge, but it is low under the negative atomic charge. Subsequently, the ion concentration is plotted in the 
opposite way in Fig. |T4](i-b) . Generating the potential gap A<f> between two ends of the system brings about 
two peculiar phenomena within the bipolar channel. At the forward bias with A<f> = IV, the electrostatic 
potential is gradually increased (Fig. 14 4i-a)), but at reverse bias with A<f> = —IV, it is sharply decreased 


(Fig. [I4[iii-a)). These two results are quite conjunctive with the ion concentration curves in the way that the 
flux is invariable along the channel axis at steady state and thus the main factor altering the potential is the 
ion distribution [27] . As plotted in Fig. [l4|^ii-b), under the forward bias, both ion species are attracted to 
the junction, so the peak value of the ion concentration is greater than that under no bias. However, under 
the reverse bias, both ion species are moved away from the junction and each one produces a small pile at 
the opposite atomic charge as presented in Fig. 14 iii-b). Accordingly, the forward bias brings about an ion 
accumulation zone at the channel junction, whereas the reverse bias creates an ion depletion zone there. 

Figure [15] shows the current-voltage curves of each ion species in the bipolar channel. While both current 
graphs are almost zero at every reverse bias, they are monotonically increased as the potential difference gets 
bigger. This result comes from the two facts that the ion-depletion zone under reverse bias terminates the 
flow inside the bipolar channel, but the ion-accumulation zone under forward bias encourages more ions to 
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Figure 14: Potential and ion concentration of a bipolar nanofluidic channel, (a) Electrostatic potential profiles; (b) 
Ionic concentration distributions of a bipolar channel along the channel axis. Three different cases include (i) no bias 
A<f> as= OV, (ii) forward bias A<f> = IV and (iii) reverse bias A<f> = —IV, where A<h = A<f># — A<E>l. We assume the 
bulk ion concentration Co to be 0.1M and the atomic charges of the channel at left and right halve, respectively, to be 
0.08e c and —0.08e c . With forward bias, both ions tend to be accumulated at the middle of the channel length, whereas 
with backward bias, both ions tend to be depleted at the junction of the channel. 
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Figure 15: Ionic current vs the applied potential difference A<I> in a bipolar nanofluidic channel. The atomic charge 
\Qk\ = 0.08e c and the bulk ion concentration Co — 0.1M are assumed. Under reverse bias, both ion species cannot 
pass through the channel. However, under forward bias, the currents of both ion species get enlarged. Especially, the 
cation current (square) increases faster according as increases. 


pass through the channel. Moreover, the positive ionic current enhances more abruptly. Another remarkable 
discovery from the current-voltage characteristics is that a higher applied voltage reduces the gradient of the 
curve, which is analytically discussed in the literature m- Therefore, by managing the external voltage, one 
can turn on and off the current through the bipolar nanochannel. 


5.2.2. Effect of atomic charge 

Figure 16 depicts the channel current in response to the applied voltage difference at \Qk\ = 0.04e c , 
\Qk\ = 0.08e c and \Qk\ = 0.12e c . Here, we set the bulk ion concentration at two reservoirs to 0.1M. In every 
case, the total current gets enlarged as Aff> becomes larger. Moreover, the rate of change of the current with 
respect to the voltage difference gets slower. The highest atomic charge amplitude, that is, \Qk\ = 0.12e c has 
the greatest amplitude of the current curve. In contrast, the lowest atomic charge amplitude does not fully 
draw both ion species at either sides under the reverse bias and thus at some negative voltage differences the 
current is nonzero. It is interesting to note that the channel current within a bipolar nanofluidic channel can 
be perfectly regulated if the atomic charges satisfy an appropriate threshold. 


5.2.3. Effect of bulk ion concentration 

We also test our bipolar channel at two different bulk ion concentration, namely, Co = 0.05M and Co = 
0.2M, and compare the total current-voltage graphs with that at Co = 0.1M as described in Fig. [TT] Every 
current-volt age curve nearly vanishes when Aff> is negative, but significantly increases when Aff> is positive. 
At higher bulk ion concentration, more ions are accumulated at the middle of the bipolar channel, so the total 
current gets bigger. Moreover, it has the maximum amplitude and gradient of the current alteration with 
respect to the voltage difference. To this end, it is expected to surmise that bulk ion concentration promotes 
the quantity of the total current through the bipolar channel. Our computational outcomes are in a good 
agreement with other numerical studies in the literature m- 

5.3. A double-well nanofluidic channel 

Finally, we consider a double-well nanofluidic channel which is named after the shape of the electrostatic 
potential curve. The electrostatic potential through the channel axis in a cylindrical channel may have several 
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Figure 16: Effects of atomic charges on the current for a bipolar nanofluidic channel. Three sets of atomic charges, 
i.e., \Qk | = 0.04e c (square), \Qk\ — 0.08e c (triangle) and \Qk\ — 0.04e c (diamond) are studied. Here, the bulk ion 
concentration Co of K + and Cl - is fixed at 0.1M. All I-V curves increase when A<f> varies from —IV to IV. Greater 
magnitude of the atomic charges results in higher channel current with forward bias, but insufficient atomic charge 
may weaken the depletion zone with reverse bias and so there is a leakage. 



Figure IT: Effect of bulk ion concentration on the channel current for a bipolar nanofluidic channel. Three sets of 
bulk ion concentrations, i.e., Co — 0.05M, Co — 0.1M and Co — 0.2M are studied. The atomic charges are given by 
\Qk \ — 0.08e c . As the bulk ion concentration gets higher, the amplitude and gradient of the current-voltage relationship 
gets maximized. 


potential wells by modifying atomic charge distribution. In fact, one of the most well-known biological channels, 
Gramicidin A channel, has a double-well transmembrane ion channel In this section, we design a 

cylindrical channel whose electrostatic potential curve has a double-well structure by varying the sign of 
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Figure 18: The schematic diagram of a 3D double-well nanofluidic channel. The channel length is divided into three 
parts. The first and last parts are negatively charged and the middle part is positively charged. The red dots indicate 
atoms with negative charges Qk = —0.12e c and the green dots indicate atoms with positive charges Qk = 0.04e c . 



Figure 19: Effect of applied voltage for a double-well nanofluidic channel, (a) Electrostatic potential profiles; (b) 
Ion concentration distributions along the double-well channel length (z— axis). Here A<f> is varied from OV (square), 
—0.2V (triangle), —0.4V (asterisk), —0.6V (diamond), —0.8V (circle) to —IV (plus sign), where A<f> is the difference 
of the applied voltage between two ends of the system. The bulk concentration is Co = 0.05M for both ions K + and 
Cl - . Two dashed vertical lines indicate the ends of the cylindrical channel. The electrostatic potential graphs shows 
two potential wells, which brings about two piles of K + ions along the channel axis. Moreover, higher applied voltage 
at the left end of the system, breaks the symmetry of the potential wells. The left well becomes weaker and the 
right well becomes stronger. Consequently, the concentration of K + ions (dashed line) at the left pile becomes lower, 
but there is little change in the concentration of K + ions at the right pile. 
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atomic charges. As illustrated in Fig.[l8j the middle section of the nanochannel is positively charged, but the 
other parts of the channel are negatively charged. 



|a<d| (V) 

Figure 20: Effect of bulk ion concentration on the channel current for a double-well nanofluidic channel. The total 
channel current versus the external voltage difference (I-V) are studied at four different bulk concentrations: Co = 0.01M 
(square), Co = 0.05M (triangle), Co — 0.1M (diamond), and Co = 0.2M (circle). With a high bulk ion concentration, 
the total current gets increased and the I-V characteristics becomes linear. 


5.3.1. Effect of applied voltage 

At first, we alter the applied voltage, but fix bulk ion concentration at Co = 0.05M and the atomic charge 
distribution as described in Fig. [18] Here, is set to be OV and is increased gradually from OV to IV. 
Figure 19 presents the electrostatic potential and ionic concentration along the channel length. On the left 
hand side of the inner channel, the electrostatic potential becomes higher, which results in moderating the left 
potential well as in Fig, [Tf^a). Subsequently, the positive ion concentration shows a dramatic change on the 
left hand side. Moreover, the small change in the potential at the right hand side of the channel corresponds 
to the small change in the concentration profile on the right. 


5.3.2. Effect of bulk ion concentration 

As in Fig.[20j the increase in the total current through the double-well channel is derived from the increase 
in the bulk ion concentration. Herein, the external voltage difference is the same (A<F = —IV). Like the 
negatively charged channel, the I-V relation becomes linear as the bulk ion concentration gets multiplied. 
These results are consistent with those observed from both numerical simulations m and experimental 
measurements m of the Gramicidin A channel. Therefore, atomic design of 3D nanofluidic channels proposed 
in the present work can be used to study biological channels, which is particularly valuable when the structure 
is not available. 


6. Concluding Remarks 

Recently the dynamics and transport of nanofluidic channels have received great attention. As a result, 
related experimental techniques and theoretical methods have been substantially promoted in the past two 
decades mm- Nanofluidic channels are utilized for a vast variety of scientific and engineering applications, 
including separation, detection, analysis and synthesis of chemicals and biomolecules. Additionally, inorganic 
nanochannels are manufactured to imitate biological channels which is of great significance in elucidating 
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ion selectivity and ion current controllability in response to an applied field in membrane channels EH EE]. 
Molecular and atomic mechanisms are the key ingredients in the design and fabrication of nanofluidic channels. 
However, atomic details are scarcely considered in nanofluidic modeling and simulation. Moreover, previous 
simulation of transport in nanofluidic channels has been rarely carried out with three-dimensional (3D) realistic 
physical geometry. Present work introduces atomistic design and simulation of 3D realistic ionic diffusive 
nanofluidic channels. 

We first proposes a variational multiscale paradigm to facilitate the microscopic atomistic description of 
ionic diffusive nanochannels, including atomic charges, and the macroscopic continuum treatment of the solvent 
and mobile ions. The interactions between the solution and the nanochannel are modeled by non-electrostatic 
interactions, which are accounted by van der Waals type of potentials. A total energy functional is utilized 
to put macroscopic and microscopic representations on an equal footing. The Euler-Lagrange variation leads 
to generalized Poisson-Nernst-Planck (PNP) equations. Unlike the hypersurface in our earlier differential 
geometry based multiscale models [94h96; , the solid-fluid interface is treated as a given profile. A domain 
characteristic function is introduced to replace the hypersurface function in our earlier formulation. 

Efficient and accurate numerical methods have been developed to solve the proposed generalized PNP 
equations for nanofluidic modeling. Both the Dirichlet-Neumann mapping and matched interface and boundary 
(MIB) methods employed to solve the PNP system in 3D material interface and charge singularity. Rigorous 
numerical validations are constructed to confirm the second-order convergence in solving the generalized PNP 
equations. 

The proposed mathematical model and numerical methods are employed for 3D realistic simulations of ionic 
diffusive nanofluidic systems. Three distinct nanofluidic channels, namely, a negatively charged nanochannel, 
a bipolar nanochannel and a double-well nanochannel, are constructed to explore the capability and impact 
of atomic charges near the channel interface on the channel fluid flow. We design a cylindrical nanofluidic 
channel of 49A in length and 10A in diameter. Several charged atoms of about 1.8 angstrom apart are equally 
located outside the channel to regulate nanofluidic patterns. For the negatively charged channel, all of the 
atoms have the negative sign; on the other hand, for the bipolar channel, half of them has the negative sign 
and the other half has the positive sign. A double-well channel has positively charged atoms at the middle 
and negatively charged atoms on the remaining part of the channel. Each end of the channel is connected 
to a reservoir of KC1 solution and both reservoirs have the same bulk ion concentration. Asymmetry in the 
applied electrostatic potentials at the ends of two reservoirs gives rise to current through these nanochannels. 
We perform numerical experiments to explore electrostatic potential, ion concentration and current through 
the channels under the influence of applied voltage, atomic charge and bulk ion concentration. 

The negatively charged channel generates a unipolar current because the negative atomic charge attracts 
counterions, but repels coions. The current within the nanochannel increases when external voltage, mag¬ 
nitude of atomic charge and/or bulk ion concentration are increased. However, the bulk ion concentration 
has a limitation in its growth because a larger bulk ion concentration shortens Debye length and thus the 
charged channel may behave like an uncharged one showing the Ohm’s law. The bipolar channel can create 
accumulation or depletion of both ions in response to the current direction. When the right end has a higher 
voltage, both ions are stored at the junction of the channel length. On the contrary, when the left end has a 
higher voltage, both ions are moved away from the junction. Applied voltage, atomic charge and bulk ion con¬ 
centration affect the amplitude and gradient of the current-volt age characteristic. At last, the special atomic 
charge distribution of the double-well channel produces the electrostatic potential profile with two potential 
wells. Increasing applied voltage at the left hand side of the system results in an obvious change in the left 
potential well and the K + concentration on the left. 

The present study concludes that the properties and quantity of the current though an ionic diffusive 
nanochannel can be effectively manipulated by carefully altering applied voltage, atomic charge and bulk ion 
concentration. Our results compare well with those of experimental measurements and theoretical analysis in 
the literature. Since the physical size of model is close to realistic transmembrane channels, the present model 
can be utilized not only for ionic diffusive nanofluidic design and simulations, but also for the prediction of 
membrane channel properties when the structure of the channel protein is not available or changed due to the 
mutation. 

Non-electrostatic interactions, are considered in our theoretical modeling but are omitted in the present 
numerical simulations to focus on atomistic design and simulation of 3D realistic ion diffusive nanofluidic 
channels. However, non-electrostatic interactions can be a vital effect in nanofluidic systems. A systematical 
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analysis of non-electrostatic interactions is under our consideration. 
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Appendices 


Table 5: Positions and charges of all atomic charges in the negatively charged channel. 


| k | x(A) y( A) z(A) Q k {e c ) || k \ x{k) y{ A) z(A) Q k (e c ) || k \ x(A) y( A) z( A) Q k (e c ) || k \ x(A) y( A) z(k) Q k (e c ) | 
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Table 6: Positions and charges of all atomic charges in the bipolar channel. 


| k | x(A) y(A) z(A) Q k (e c ) || k \ x(A) y( A) 2(A) Q k {e c ) || k \ x{A) y( A) z( A) Q k (e c ) || k \ x(A) y{A) z{ A) Q k (e c ) \ 
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